% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%latest date 30 OCT 2005
clear;
data = load('c:\manoj\projects\eej\eej_par_all.txt');
%
%data1 = load('c:\manoj\projects\eej\eejdata.txt'); %<ettmnhyb 2000-2001; 
%data1 = load('c:\manoj\projects\eej\eejdata_new.txt');%<ettmnhyb2000-2002
%and corrected tirmnabg

%----31 Jan 2006 Added HUA - FUQ test 
%data1 = load('c:\manoj\projects\eej\eejdata31jan.txt');
%mjd,y,m,d,j-1,ettmnhyb(pp),ettmnabg(pp),tirmnhyb(pp), tirmnabg(pp),huamnfuq(pp),kp);

%-------------------------------------------------------

%---1Feb added MBO-GUI
data1 = load('c:\manoj\temp\eej_obs_data.txt');

%15 Jan 2006 Some notes
% Basically tried to veryfy the results. Checked the time stamps ->
% corrected the IST-Solar Time problems.. Checked with many LT bin
% combinations
% someof the resluts are mentioned in correlation.doc 
% LT 9-15 and bins 10 degree with the necessary corrections brings the most
% aestetic (?) results

%-------- correlate the CHAMP derived EEJ strength and ground delta H
% without considering the lT effect
%--------------------------------------------------------------------

% nloop = 1
% for ki = -33:10:133
% L = data(:,2) > ki & data(:,2) < ki+10;
% IndiaE = data(L,:);
% for i = 1:length(IndiaE(:,1)),
% % L = data1(:,1) < IndiaE(i,1);
% % ind(i) = sum(L);
% ind(i) = findnearest(IndiaE(i,1),data1(:,1));
% end;
% obsE = data1(ind,:);
% L = obsE(:,7) > 500;
% obsE(L,:) = [];
% IndiaE(L,:) = [];
% 
% cof = corrcoef(IndiaE(:,9),obsE(:,7));
% ceff(nloop) = cof(2,1);
% nloop = nloop +1;
% clear ind;
% end;

%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------

%Control Center

sector = 'afr1';%options 'ind' 'amr' 'afr' 'afr1' 'pac'
obsdata = 7;  %ETT-HYB = 1,ETT-ABG = 2;TIR-HYB = 3;TIR-ABG = 4; HUA-FUQ 5; MBO-GUI = 6; AAE-QSB 7; GUA-CBI 8;
satdata = 2; % %Pear Current Density = 1; av cur distan 2;
obslt = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
latbin = 10; %5,10,20
plotoptions = 'b.-';
%<------- one problem  assumes that TIR is the center (77.8)

switch sector,
    case 'ind',
        obslon = 82.5; %This is the longitude of ALLahabad. The difference beteen allahabad and actual longitude of the observatory is taken care
        switch latbin,
    case 5,
        x1data = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
        incrlat = 5;
        st = 15.3; 
        en = 135.3;
    case 10,
        x1data = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
        incrlat = 10;
        st = 12.8;
        en = 132.8;
    case 20,
        x1data = [-60 -40 -20 0 20 40 60];
        incrlat = 20;
        st = 7.8;
        en = 147.8;
    otherwise,
        display('Error');
        break;
        return;
    end;
    
    case 'amr',
        obslon = 285;
    switch latbin,
    case 5,
        x1data = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
        incrlat = 5;
        st = -137.5; 
        en = -17.5;
    case 10,
        x1data = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
        incrlat = 10;
        st = -139.5;
        en = -19.5;
    case 20,
        x1data = [-60 -40 -20 0 20 40 60];
        incrlat = 20;
        st = 7.8;
        en = 147.8;
    otherwise,
        display('Error');
        break;
        return;
        end;
 
    case 'afr',
        obslon = 343;
    switch latbin,
%     case 5,
%         x1data = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
%         incrlat = 5;
%         st = -137.5; 
%         en = -17.5;
    case 10,
        x1data = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
        incrlat = 10;
        st = -81.5;
        en = 47.5;
%     case 20,
%         x1data = [-60 -40 -20 0 20 40 60];
%         incrlat = 20;
%         st = 7.8;
%         en = 147.8;
    otherwise,
        display('Error');
        break;
        return;
        end;        

 case 'afr1',
        obslon = 37.5;
    switch latbin,
%     case 5,
%         x1data = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
%         incrlat = 5;
%         st = -137.5; 
%         en = -17.5;
    case 10,
        x1data = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
        incrlat = 10;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
%     case 20,
%         x1data = [-60 -40 -20 0 20 40 60];
%         incrlat = 20;
%         st = 7.8;
%         en = 147.8;
    otherwise,
        display('Error');
        break;
        return;
        end;            

case 'pac',
        obslon = 143;
    switch latbin,
%     case 5,
%         x1data = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
%         incrlat = 5;
%         st = -137.5; 
%         en = -17.5;
    case 10,
        x1data = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
        incrlat = 10;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
%     case 20,
%         x1data = [-60 -40 -20 0 20 40 60];
%         incrlat = 20;
%         st = 7.8;
%         en = 147.8;
    otherwise,
        display('Error');
        break;
        return;
        end;                    
        
end;


switch obsdata,
    case 1,
        load c:\manoj\m\eej\17jan06 ETTMNHYB;
        k = nanmean(ETTMNHYB);
        k=circshift(k,[0,5])'; % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 7;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
    case 2,
       load c:\manoj\m\eej\17jan06 ETTMNABG;
        k = nanmean(ETTMNABG');
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 8;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
    case 3,
       load c:\manoj\m\eej\17jan06 TIRMNHYB;
        k = nanmean(TIRMNHYB');
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 9;
        data1(:,1) = data1(:,1) + ((77.85-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78.5 (TIR) and 82.5 (Allahabad - IST)
        xtime = 1:24;
    case 4,
        load c:\manoj\m\eej\17jan06 TIRMNABG;
        k = nanmean(TIRMNABG');
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 10; 
        data1(:,1) = data1(:,1) + ((75-82.5)/15)/24; % This assumes 75 (between TIR and ABG)
        xtime = 1:24;
    case 5,
        load c:\manoj\projects\eej\HUAMBOAAE EEJindex;
        k = nanmean(EEJindex,2);
        k = circshift(k,19);
        xtime = 0.5:1:23.5;       
        obs_index = 11;
    case 6,
        load c:\manoj\projects\eej\HUAMBOAAE MBOMNGUI;
        k = nanmean(MBOMNGUI)';
        k = circshift(k,23); % 23 ~ = 343/15
        xtime = 0.5:1:23.5;       
        obs_index = 12;
    case 7,
        load c:\manoj\projects\eej\HUAMBOAAE AAEMNQSB;
        k = nanmean(AAEMNQSB)';
        k = circshift(k,2); % 2 ~ = 37.5/15
        xtime = 1:24;       
        obs_index = 13;        
    case 8,
        load c:\manoj\projects\eej\HUAMBOAAEGUA GUAMNCBI;
        k = nanmean(GUAMNCBI)';
        k = circshift(k,9); % 9 ~ = 144/15 (the remaining would be taken care by shifting time axis to 1:24 (instead 0.5:1:23.5)
        xtime = 1:24;       
        obs_index = 14;                
    otherwise,
        display('Error');
        break;
        return;
end;

switch satdata,
    case 1,
        eej_index = 9;
    case 2,
        eej_index = 10;
    otherwise,
        display('Error');
        break;
        return;
end;        



[p,s] = polyfit(xtime',k,9);
%----------------------------

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;

for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < max(data1(:,1)); %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data1(ind,1);
    chpt(i) = IndiaE(i,1);
    
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 , %The last fliter take care the nodata (9999) points
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,lt1); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,lt2); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    [cof,err] = corrcoef(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)');  % correlation & Correction for observatory data 
%   [cof,err] = corrcoef(IndiaE(:,eej_index),obsE(:,obs_index));
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;
    
    
    %----------- PLOT SCATTER
% subplot(2,4,nloop-1);
% plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)','rs','MarkerFaceColor','r');
% axis([0 0.3 0 200]);
% xlabel([sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ' 'Coef ' sprintf('%4.2f',cof(2))]);
% 
%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt;
end;

x1data(isnan(x1data)) = [];
fprintf('\n');
hold on;
errorbar(x1data,ceff,cerr/2,plotoptions);
axis([-80 80 -.6 1.4]);
grid on;

% %Notes
% %16.1.2006 Checked KP of all the selected events all are < 2.0 as expected
% %since the satellite data are screened for kp
% % lt2= ltt2;
% % lt2(ind2) = [];
% % L = lt2 >= 10 & lt2 <= 13;
% subplot(421);
% plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)','rs','MarkerFaceColor','r');
% %hold on;plot(IndiaE(L,eej_index),obsE(L,obs_index).*(f_ch(L)./f_ob(L))','bs','MarkerFaceColor','b');
% xlabel('Pear Current Density (CHAMP)');
% ylabel('\delta H nT');
% axis([0 0.5 0 200]);
% %legend(['9-15: Coef ' sprintf('%4.2f',cof(2))],['10-3: Coef ' sprintf('%4.2f',cof1(2))]);
% legend(['10-3: Coef ' sprintf('%4.2f',cof(2))]);
% % title(['ETT HYB  ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
% % 
